% By Martin M. Andreasen, April 22 2010
% This function simulates the model when solved up to third order. The purning scheme is used.
% These codes are faster for a third order approximation because we use the kronecker
% product to eleminate loops
function [y_sim,xf_sim,xs_sim,xrd_sim] = simulate_3rd_kron_v2(gx,gxx,gxxx,hx,hxx,hxxx,gss,gssx,gsss,hss,hssx,hsss,eta,shocks,orderApp)

% Some indices
ny      = size(gx,1);        %Number of output variables
nx      = size(hx,1);        %Number of state variables
ne      = size(eta,2);       %Number of structural shocks
num_sim = size(shocks,2);    %Number of simulations

% Allocating memory
y_sim   = zeros(ny,num_sim);
xf_sim  = zeros(nx,num_sim);
xs_sim  = zeros(nx,num_sim);
xrd_sim = zeros(nx,num_sim);

% The simulation
xf    = zeros(nx,1);    %The first order terms
xs    = zeros(nx,1);    %The second order terms
xrd   = zeros(nx,1);    %The third order terms

% Defining matrices
HHxxtil  = 1/2*reshape(hxx,nx,nx^2);
HHxxxtil = 1/6*reshape(hxxx,nx,nx^3);
GGxxtil  = 1/2*reshape(gxx,ny,nx^2);
GGxxxtil = 1/6*reshape(gxxx,ny,nx^3);

if orderApp == 1
    for t=1:num_sim
        % The controls
        xf_sim(:,t) = xf;
        y_sim(:,t) = gx*xf;
        
        % The state variables
        xf_p  = hx*xf + eta*shocks(:,t);
        
        % Updating xf
        xf  = xf_p;
    end
elseif orderApp == 2
    for t=1:num_sim
        % The controls
        xf_sim(:,t) = xf;
        xs_sim(:,t) = xs;
        AA          = kron(xf,xf);
        y_sim(:,t)  = gx*(xf+xs) + GGxxtil*AA + 1/2*gss;
        
        % The state variables
        xf_p  = hx*xf + eta*shocks(:,t);
        xs_p  = hx*xs + HHxxtil*AA + 1/2*hss;
        
        % Updating xf, xs
        xf  = xf_p;
        xs  = xs_p;
    end
elseif orderApp == 3
    for t=1:num_sim
        % The controls
        xf_sim(:,t)  = xf;
        xs_sim(:,t)  = xs;
        xrd_sim(:,t) = xrd;
        AA    = kron(xf,xf);
        BB    = kron(xf,xs);
        y_sim(:,t) = gx*(xf+xs+xrd) + GGxxtil*(AA+2*BB) + 1/2*gss + ...
            GGxxxtil*kron(xf,AA) + 3/6*gssx*xf + 1/6*gsss;
        
        % The state variables
        xf_p  = hx*xf + eta*shocks(:,t);
        xs_p  = hx*xs + HHxxtil*AA + 1/2*hss;
        xrd_p = hx*xrd + 2*HHxxtil*BB + HHxxxtil*kron(xf,AA)+ 3/6*hssx*xf + 1/6*hsss;
        
        % Updating xf, xs, xrd
        xf  = xf_p;
        xs  = xs_p;
        xrd = xrd_p;
    end
end
